% 求解简单的一维弹性力学问题
% Gitee Repo

clc
clear

% 长度
L = 1;

% 杨氏模量
E = 10;

dx = 0.1;
x = (0:dx:L)'
n = size(x,1);

% 各节点位移
u = zeros(n,1);

% 各节点所受外力(体力)
b = zeros(n,1);
b = sin(pi/L*x);
%b(n)=1;

% 计算刚度矩阵
% 计算方法：
% 应力平衡方程 \dv{\sigma}{x} + f = 0
% 应变几何方程 \varepsilon = \dv{u}{x}
% 本构关系 \sigma = E \varepsilon
% 因此, E \dv[2]{u}{x} + f = 0
% 离散化，E \frac{u(n+1) - 2*u(n) +u(n-1)}{(dx)^2}+f=0
% 即 u(n+1) - 2*u(n) +u(n-1) =-f(dx)^2/E
% 边界条件：
% 第一个节点位移 u(1) = 0
% 最后一个节点位移的导数 du(n)/dx=b(n)
% 为了满足边界条件，需要特殊处理第一个和最后一个节点
K = spalloc(n,n,3*n);
K(1,1)=1;
K(n,n)=-1;
K(n,n-1)=1;
for i = 2:n-1
  K(i,i) = -2;
  K(i,i-1) = 1;
  K(i,i+1) = 1;
end

_b = -b*(dx)^2/E;
% 由于假定左端固定，因此不应对该节点重复施加外力
_b(1) = 0;
% 计算位移 F=Ku, u=inv(K)*F
u = gmres(K,_b,[],[],100);

% 计算变形
ux = zeros(n,1);
ux(1:n-1) = u(2:n)-u(1:n-1);
ux(n) = u(n)-u(n-1);
ux(1) = u(2)-u(1);

figure();
hold on
plot(x,u);
plot(x,ux);
xlabel('x')
ylabel('u')
legend('位移','变形')

figure()
hold on
axis equal
scatter(x, zeros(n,1),'k')
scatter(x+u, zeros(n,1))
line([x(1) x(n)+u(n)],[0 0])
xlabel('x')

